Applicable Algebra in Engineering, Communication and Computing manuscript No. 

(will be inserted by the editor) 



o 
o 

(N 
Q 

Q 



c/3 



> 
00 



o 



X 



Michel Fliess • Cedric Join 
Mamadou Mboup 



Algebraic change-point detection 



Received: date / Revised: date 



. . Abstract Elementary techniques from operational caleulus, differential al- 

' gebra, and noncommutative algebra lead to a new approach for change-point 

detection, which is an important field of investigation in various areas of ap- 
plied sciences and engineering. Several successful numerical experiments are 



Q I presented. 
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1 Introduction 



Let / : M ^ R be a piecewise smooth function with discontinuities at 
I ti,t2, ■ ■ ■ ■ Its pointwise derivative /'•^•' which exists and is continuous ex- 

cept at ti,t2, ■ ■ ■ , and its distribution derivative /' in Schwartz's sense are, 
as well known, related by 



f{t) = f^'Ht) + if{tl+) - f{tl-))S{t-h) + {f{t2+) - f{h~))S{t-t2) + .. . 
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where 

— /(t+) = lim^, fit), f{T-) = limtt, fit), 

— S is the Dirac delta function. 

A huge Hteratur^B has been devoted to the detection of ti , t2 , • ■ • , which is a 
crucial question in signal processing, in diagnosis, and in many other fields of 
engineering and applied sciences, where those discontinuities are often called 
change-points or abrupt changesiQ Difficulties are stemming from 

— corrupting noises which are blurring the discontinuities, 

— the combined need of 

— fast online calculations, 

— a feasible software and/or hardware implementation. 

Most of the existing literature is based on statistical tools (see, for instance, 
[UmiSllH] and the references therein). 

The origin of our algebraic viewpoint lies in the references p!8lfT9] which 
are devoted to the parametric identification of linear systems in automatic 
control0 We employ elementary techniques stemming from operational cal- 
culusQ differential algebra and noncommutative algebra. We are replacing 
Eq. ^ by its operational analogue which is easier to handle. Restricting 
ourselves to solutions of operational linear differential equations with ratio- 
nal coefficients lead to noncommutative rings of linear differential operators. 
By representing a change-point by a delay operator, i.e., an operational ex- 
ponential. Sect. [2] concludes with the identifiability of change-points, i.e., 
the possibility of expressing them via measured dataQ Higher order change- 
points, i.e., discontinuities of derivatives of various orders are briefly discussed 
in Sect.O Sect. [4] presents several successful numerical experiments!! which 

— exhibit good robustness properties with respect to several types of addi- 
tive and multiplicative corrupting noises; 

— indicate that our approach is still valid outside of its full mathematical 
justification0 

Preliminary results may be found in [17] and p!6l[27] . 

^ See the excellent account due to Basseville and Nikiforov [l] for more details. 
^ The most popular terminology in French is ruptures. 

^ Ch ange-point detection has also been studied in [5] via tools stemming from 
[181119] . but in a quite different manner when compared to us. 

* Mikusinski's foundation [291130] of operational calculus, which is not based on 
the usual Laplace transform, is a better choice for the connection with the other 
algebraic tools. Mikusinski's work, which is a superb example of algebratc analysis, 
is to o much neglected in spite of some advertisements like the nice book by Yosida 

m- 

^ In the context of constant linear control systems with delays, which bears some 
simil arity w ith what is done here , the id entification of delays has also been tackled 
in [3l l3T]|36] via techniques from [T8][T9| . 

® Let us emphasize that our techniques have already been applied in some con- 
crete ca se-stud ies, where the sig nals to be processed are stemming from either 
biology [^157] or finance [IMI3- 

^ It goes without saying that this Section, which is mainly descriptive, is not 
intended to be fully rigorous. 
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2 Algebra via operational calculus 

2.1 Differential equations 

Take a commutative field kg of characteristic zero. The field fco(s) of rational 
functions over ko in the indeterminate s is obviously a differential field with 
respect to the derivation ^ and its subfield of constants is /co (cf. [6[l34j). 
Write fco(s)[^] the noncommutative ring of linear differential operators of 
the form 

d 



d" 



ds 



Qa e ko{s) (2) 



We know that fco(s)[^] is a left and right principal ideal ring (cf. [28ll34pl ). 
Any signal x is assumec0 here to be operationally holonomic, i.e., to satisfy 
a linear differential equation with coefficients in ko{s): there exists a linear 
differential operator zu G fco(s)[^], 7^ 0, such that wx — 0. 

Remark 1 Let us explain briefly this assumption. We consider only holo- 
nomic time functions z{t)^ i.e., time functions which satisfy linear differential 
equations with polynomial coefficients: 

The corresponding operational linear differential equation reads (cf. '39^) 




where / G C[s] depends on the initial conditions. A homogeneous linear 
differential equation is obtained by differentiating both sides of the previous 
equation enough times with respect to s. 

Let K be the algebraic closure of fco(s): K is again a differential field with 
respect to ^ and its subfield of constants is the algebraic closure k^ of /cq- It 
is known that x belongs to a Picard-Vessiot extension of K (cf. [51l34j). 

Remark 2 Holonomic functions play an important role in many parts of 
mathematics like, for instance, combinatorics (see, e.g., |10j). 

* Note that [M] is not employing, contrarily to j28j , the usual terminology of ring 
and module theory. 
^ See also [l2l[lfl|25] . 



4 



Fliess, Join & Mboup 



2.2 Annihilators 

Consider now the left A;o('S)[^]-niodule M spanned by a finite set {xjt G 1} 
of such signals. Any Xi, is a torsion element (cf. [53]) and therefore is 
a torsion moduler^l The annihilator Aj of {xji G /} is the set of linear 
differential operators rui G ^o('S)[^] such that, \/ l £ I, wix^ — 0. It is 
a left ideal of fco(s)[^] and it is therefore generated by a single element 
A G fc(s)[^], Z\ 7^ 0, which is called a minimal annihilator of {xji G /}. It is 
obvious that A is annihilating any element belonging to the fco"Vector space 
spanj,^(xt|i G /). The next result is straightforward: 

Lemma 1 Let Ai, A2, Ai ^ A2, be two minimal annihilators. There exists 
p d ko{s), p ^ 0, such that Ai = pA2. 

We will say that the minimal annihilator is unique up to left multiplications 
by nonzero rational functions. 

A rational function ^, p,q G ko[s], q ^ 0, is said to be proper (resp. 
strictly proper) if, and only if, the d°p < d°q (resp. d°p < d°q). A differential 
operator ^ is said to be proper (resp. strictly proper) if, and only if, any ^q, 
is proper (resp. strictly proper). The next result is an obvious corollary of 
Lemma [TJ 

Corollary 1 It is possible to choose an annihilator, which is minimal or not, 
in such a way that it is proper (resp. strictly proper). 

A rational function ^, p,q ko[s], q ^ 0, is said to be in a finite integral 
form (resp. strictly finite integral form) if, and only if, it belongs to fco[-j] 
(resp. A differential operator Q is said to be in a finite integral 

form (resp. strictly finite integral form) if, and only if, any pa is in a finite 
integral form (resp. strictly finite integral form). Consider a common multiple 
m G k[s] of the denominators of the QaS. The operator s~^m-uj is in a 
(strictly) finite integral form for large enough values of the integer > 0. 

Corollary 2 It is possible to choose an annihilator, which is minimal or not, 
in such a way that it is in a finite integral form (resp. strictly finite integral 
form). 

2.3 Delay operators 

Let k/ko be a transcendental field extension. The field k{s) of rational func- 
tions over k in the indeterminate s is again a differential field with respect 
to and its subfield of constants is k. The noncommutative ring k{s)[j^] 
of linear differential operators is defined as in Sect. 12.11 Pick up an element 
tr G k, called delay, which is transcendental over ko. Write the delay operator 
with its classic exponential notation e"'''" (cf. [33]), as it satisfies the differ- 
ential equation -I- tr) e"*''" = 0. According to Sect. 12.21 the differential 
operator + G A;(s)[^] is a minimal annihilator of e"*''*. 



Such a module is called a differential module in [34j . 
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2.4 Identifiabihty of the delay 
2.4-1 Main result 

Let vui, wi G ^o(s)[^] be minimal annihilators of two signals xi, X2, x\X2 ^ 
0. Introduce the quantity 

X = xi + xae-*"" (3) 

Multiplying on the left both sides of Eq. ([3]) by wx yields vd\X = vD\X2e~'*"^^ . 
Thus 

tuiXe*'''' = vjxx^ 

and 

ro^TUiXe*'-'* = (4) 

where G ^o(s)[^] is a minimal annihilator of tnia;2. The next proposition 
follows at once: 

Proposition 1 Eq. ([4]) is equivalent to 

^ tj: (tt^X) = 0, 7r^Gfco(s) 

finite 

where at least one tTu is not equal to 0. 
Write ko{s){X) the differential overfield of fco(s) generated by X. 
Corollary 3 tr in Eq. ([3|) is algebraic over the differential field ko{s){X). 

Remark 3 Assume that the quantity X is measured, i.e., there exists a sensor 
which gives at each time instant its numerical value in the time domain. Then, 
according to the terminology in J15], Corollary [3] may be rephrased by saying 
that tr is algebraically identifiable. 

2.4-2 First example 

Set Xi — 'Y^^i=o gJ^^i ; * = 1:2, 7i,. e kfj, where Ni is a known non-negative 

integer l"1 Then s^'^^ is a minimal annihilator of Xi. It follows at once 

that Proposition |1] and Corollary [3] apply to this case. 

Straightforward calculations demonstrate that tr is the unique solution 
of an equation of the form 

p(^^+tr) X^Q (5) 

where p e fco(s)[f ], 1 < g < Af2. 

The coefficients 71^^ are not necessarily known. 
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2.4-3 Second example 

Assume that X2 — N^+i i 72 G ^o, in Sect. 12.4.21 Multiply both sides of Eq. 
m by if,+U)s^'+^ yields 

as as 

Eq. ^ becomes 

M-^+-tr)s'''+'X = (6) 

where tti G A:o(s)[^] is a minimal annihilator of ( + tr)s'^^~^^xi. 
Proposition 2 t,. satisfies an algebraic Equation ^ of degree 1. 

Remark 4 If X is measured as in Remark|31 then, according to the terminol- 
ogy in J18','19] , Proposition [2] may be rephrased by saying that tr is linearly 
identifiable. 

2.4-4 Third example 

Assume in Eq. ([3]) that X2 = f G kois), a. b £ ko[s], (a, 6) = 1, is a known 
rational function, i.e., 

X^x, + ^e-*-^ (7) 

Multiplying both sides by ijs+U)^ yields {-^+tr)^X = {d.+tr)^xi. Since 
tr is constant, there exists an annihilator tt G fco(s)[^] of + tr)^xi, i.e., 

Proposition 3 tj. in Eq. ([7]) satisfies an algebraic Equation ([5]) of degree 1. 

Remark 5 If X is measured as in Remark [31 then, according to Remark 21 
Proposition [3l may be rephrased by saying that tr is linearly identifiable. 

3 Higher order change-points 

Take again as in the introduction a piecewise smooth function /, which is 
now assumed to be C", n > 0, i.e., / and its pointwise derivatives up to 
order n are continuous. We might be interested in the discontinuities of its 
{n + iy^ order pointwise derivative, which are called change-points, or abrupt 
changes, of order n + I. 
By replacing Eq. ^ by 

it is straightforward to extend all the results of Sect. 12.41 to higher order 
change-points. 
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4 Some numerical experiments 

4.1 General principles 

From now on fco is a subfield of R, Q for instance. We utilize the calculations 
of Sect. [US] like follows: 

— Multiplying both sides of Eq. ^ by s^^ , where > is large enough, 
yields 

S-^TTl (^^+t,^5^^ + lX = (9) 

where s~^Tri{j^ + tr)s^^'^^ is a strictly integral operator. 

— Going back to the time domain is achieved via the classic rules of oper- 
ational calculus P^f5Ul[55] . where j-^ corresponds to the multiplication 
by i-tr^ 

— xi = X]i/i=o T^TTT and X2 = gN2+i correspond in the time domain to the 

polynomial functions E^^Lo -^Trf- ^nd ^^2^^. 

— Those time functions are assumed to approximate on a "short" time in- 
terval the signal where change-points have to be detected. 

— Consider the numerical value v taken by the time analogue of the left 
side of Eq. ([9]) when the value given to U is the middle of a given "short" 
time window. If v is "close" to 0, then we say that the middle of the time 
window is a change-point. 

— This time window is sliding in order to capture the various change-points, 
which are assumed to be not too "close", i.e., the distance between two 
consecutive change-points is larger than the time window. 

— The corrupting noises are attenuated by the iterated time integrals which 
corresponds in the time domain to the negative power of s in the left side 
of Eq. dUH 



4.2 ExamplefO 

The following academic examples are investigated: 

— piecewise constant and polynomial real-valued functions, 

— a real-valued sinusoid plus a piecewise constant real-valued function. 



Noises in [Tl] are viewed, via nonstandard analysis, as quickly fluctuating phe- 
nomena (see also |23) for an introductory presentation). The noises are attenuated 
by the iterated time integrals, which are simple examples of low-pass filters (we 
may also choose, according to Lemma (2] more involved low-pass filters (see, e.g., 
[3))- No statistical tools are required and we are by no means restricted to Gaus- 
sian white noises, like too often in the engineering studies. Moreover the corrupting 
noises need not to be additive. They might also be multiplicative. 

Interested readers may ask C. Join for the corresponding computer programs 
(Cedric . Join@cran.uhp-nancy . f r). 
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The robustness with respect to corrupting noises, which is reported in Table 
1, is tested thanks to several noises, of various powersP^ which are of the 
following types: 

1. additive, zero mean, and either normal, uniform or PcrlinF^ 

2. multiplicative, of mean 1, and uniform. 

We finally note that 

— piecewise polynomial functions were difficult to analyze even via recent 
techniques hke wavelets (see, e.g., [9l[24]): 

— we do not need any a priori knowledge of the upper bound of the number 
of change-points (see, e.g., [2T|[35]): 

— we are not limited to a given type of noises and we are able to handle 
multiplicative noises as well (see, e.g., [20l[22ll38] ) : 

— the results remain satisfactory even with a very high noise level (see Fig- 
ures [1] and [2]) . 

Remark 6 The so-called Perlin noises, which are not familiar in signal pro- 
cessing and in automatic control, contain components which are obviously 
not quickly fluctuating. It is all the more remarkable that our computer sim- 
ulations are still good, in spite of the fact that this example goes beyond the 
theoretical justifications provided in Sect. 14.11 




50 100 150 200 250 300 350 400 450 500 50 100 150 200 250 300 350 400 450 500 



(a) Noise-free signal (- -), signal (-) (b) Change-point detection - Exact 
Fig. 1 Piecewise constant signal - Normal additive noise - SNR: db 



We are utilizing the notion of signal-to-noise ratio, or SNR, which is familiar 
in signal processing (see Wikipedia, for instance). 
^■^ Perlin's noises [32] are quite popular in computer graphics. 
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50 100 150 200 250 300 350 400 450 500 50 100 150 200 250 300 350 400 450 500 



(a) Noise-free signal (- -), signal (-) (b) Change-point detection - Exact (+) 
Fig. 2 Piecewise constant signal - Normal additive noise - SNR: -6 db 




500 1000 1500 2000 2500 3000 3500 4000 500 1000 1500 2000 2500 3000 3500 4000 

(a) Noise-free signal (- -), signal (~) (b) Change-point detection - Exact (+) 
Fig. 3 Piecewise polynomial signal - Normal Additive noise - SNR: 25 db 




500 1000 1500 2000 2500 3000 3500 4000 



500 1000 1500 2000 2500 3000 3500 4000 



(a) Noise-free signal (- -), signal (-) (b) Change-point detection - Exact 
Fig. 4 Piecewise polynomial signal - Uniform additive noise - SNR: 25 db 
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500 1000 1500 2000 2500 3000 3500 4000 




500 1000 1500 2000 2500 3000 3500 4000 



(a) Noisy free signal (- -), signal (-) (b) Change-point detection - Exact (+) 
Fig. 5 Piecewise polynomial signal - Additive Perlin noise - SNR: 20 db 




500 1000 1500 2000 2500 3000 3500 4000 500 1000 1500 2000 2500 3000 3500 4000 



(a) Noise-free signal (- -), signal (-) (b) Distribution of change-point detec- 
tion, exact places 

Fig. 6 Piecewise polynomial signal - Uniform multiphcative noise - SNR: 20 db 
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500 1000 1500 2000 2500 3000 3500 4000 4500 5000 

(a) Noise-free signal (- -), signal (-) (b) Change-point detection - Exact 
Fig. 7 Sinusoidal signal - Normal additive noise - SNR: 25 db 
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500 100015002000250030003500 4000 4500 5000 500 100015002000250030003500 4000 4500 5000 

(a) Noise-free signal (- -), signal (-) (b) Change-point detection - Exact 
Fig. 8 Sinusoidal signal - Normal additive noise - SNR: 20 db 
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500 1000 1500 2000 2500 3000 3500 4000 4500 5000 



500 1000 1500 2000 2500 3000 3500 4000 4500 5000 



(a) Noise-free signal (- -), signal (-) (b) Change-point detection - Exact (-h) 
Fig. 9 Sinusoidal signal ~ Additive Perlin noise - SNR: 10 db 
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Table 1 Summary of simulation results • -|-: additive noise • x: multiplicative 
noise 
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